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Abstract 



We perform a two-dimensional molecular-dynamics study of a model for 
sheared bidisperse granular systems under conditions of simple shear and 
Poiseuille flow. We propose a mechanism for particle-size segregation based on 
the observation that segregation occurs if the viscous length scale introduced 
by a liquid in the system is smaller than of the order of the particle size. We 
show that the ratio of shear rate to viscosity must be small if one wants to find 
size segregation. In this case the particles in the system arrange themselves 
in bands of big and small particles oriented along the direction of the flow. 
Similarly, in Poiseuille flow we find the formation of particle bands. Here, 
in addition, the variety of time scales in the flow leads to an aggregation of 
particles in the zones of low shear rate and can suppress size segregation in 
these regions. The results have been verified against simulations using a full 
Navier-Stokes description for the liquid. 

PACS numbers: 83.20.Hn, 47.55.Kf, 02.70.Ns 
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I. INTRODUCTION 



If a mixture of particles and liquid is sheared, either in a continuous fashion or by periodic 
excitation of the system, a host of structural rearrangements in the mixture are known to 
occur. Most of the systems studied experimentally involve mono-disperse suspensions . 
For concentrated suspensions and oscillatory shearing the formation of layers in the iso- 
velocity planes and an additional arrangement in lines of particles perpendicular to the 
velocity field |J has been observed. In addition, simulations of Taylor- Couette flow in 
the viscous regime confirm the organization of the particles in the suspension in layers 
oriented along the iso- velocity planes ||. Numerically, cluster formation has been found 
in these simulations within the formed planes. Although for bidisperse suspensions similar 
organization patterns have been found [RJ and size-distribution induced instabilities have 
been observed in sedimentation experiments [[J, simulations here so far have not found clear 
indications of size segregation under shear. 

In this paper, we want to investigate particle size segregation of a granular medium 
under shear flow conditions in two dimensions in the presence of a liquid. The mechanism 
for size segregation in this model is basically of geometric nature and may generalize to three 
dimensions. We concentrate on pipe and laminar shear flows described by local shear rates 
7 = (d/dy)v x , where we have in mind (see Fig. |T]) that the flow is translationally invariant 
in x direction and that its speed varies in y direction. 

II. THE MODEL 

For the sake of simplicity we assume the particle size distribution to be bidisperse, i.e., 
our system consists of a number n s of small and a number of big disk shaped particles 
with radii r s and r?,, respectively. These particles are initially distributed randomly over a 
two dimensional plane with a specified area fraction c. The overall area fraction c is the 
sum of the area fractions of big and small particles, c = c s + q,. Under viscous conditions 
and when two particles are nearly in contact, lubrication effects will introduce a strong 
velocity dependent damping that inhibits particle collisions. In our model, we allow for 
small particle overlaps, but between overlapping particles strong repulsive forces and and 
a velocity proportional damping are introduced. I.e, in normal direction acts a damped 
harmonic force f n , 

fn = -fj£ + 2/IKV n , (1) 

where /i is the reduced mass of the pair, £ = |xi — x 2 | — (n + r 2 ) is the virtual overlap of 
the two particles located at xx and x 2 (Fig. 0), v n is the normal component of the relative 
velocity, and k parameterizes a velocity proportional damping. The form of Eq. (|I]) is 
such that the restitution coefficient e, the ratio of the normal velocities after and before 
the impact, is the same for all binary collisions, independent of mass and velocity. The 
restitution e is related to k by k e l/\Jl + (it / log(e)) 2 and is chosen to be 0.8 in this 
work. Equation ([!]) is in dimensionless form, the chosen length scale is the average radius 
f = (n^b + n s r s )/(nt, + n s ), the mass unit is the mass of a particle of average radius, and 
the time unit is such that the "spring constant" of a pair contact becomes unity and thus 
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does not appear in Eq. (]T|). More precisely, in these units an elastic contact with k = 
has a duration of it. We neglect tangential forces in our model and consequently we do not 
consider particle rotation. The equations of motion are integrated using a fourth order Gear 
predictor-corrector algorithm, with the time step chosen to be 0.15, which still guarantees 
good energy conservation in the elastic case. 

To prepare the initial configuration, we start from a random configuration of disks and 
simulate with periodic boundary conditions until time t = 20 using a very low restitution 
coefficient of 0.1 to remove particle overlaps and excess energy in the system efficiently. 
Afterwards the setup is continued with e = 1 until t = 100, corresponding to elastic collisions. 
The average energy per particle in the system during this stage is quite low and allows only 
marginal overlaps. Thus, the resulting configuration is close to that of a hard disk system 
in thermal equilibrium. 

The liquid motion is approximated by an invariant velocity field in x direction, 

u(y) = yje x , (2) 

where y is counted from the center of the simulation cell (cf. Fig |1]). We have tested the 
implicit assumption that the velocity profile is stationary for shear flow, i.e., that it is not 
modified by the presence of the particles in the fluid @. It differs only little from the 
theoretical constant viscosity profile in the case of Poiseuille flow as shown later in this 
article. 

The liquid exerts an additional drag on each particle i which is added to the interaction 
force ([I]) between particles. The force is here assumed to be the Stokes drag force on a 
sphere, 

U,i = -67^77 [vj - u(yi)] , (3) 

where is the velocity of particle % and yi its y coordinate. 

We then start the simulation by choosing the initial particle velocities to be equal to 
the local liquid-background velocity. The boundary conditions for the particles are Lees- 
Edwards boundary conditions (as displayed in Fig. |], cf. also [0]) with shear rate 7. 

III. RESULTS OF THE SHEAR FLOW CASE 

In Figs. |3| and |] we show two typical simulation sequences whose physical parameters 
differ only in the value for the background viscosity 77. The ratio of the particle radii in our 
bidisperse system is r\,/r s = 4, the ratio of the number density of big and small particles 
is rib/n s = 0.05 and the total area fraction of particles is rather high, c = 0.6. In the long 
time limit one sees very different structural reordering of the systems emerging. In the 
first sequence (Fig. |3|) with low viscosity, the particle arrangement is more or less random, 
whereas in the second case (Fig. |4]) one observes a very clear separation into alternating 
zones parallel to the flow direction which contain alternatingly big and small particles, 
respectively. 

The basis to understand this segregation phenomenon lies in an analysis of the length 
scales that are present in the system. Apart from the system size and the two different 
particle radii an additional viscous length exists in the problem. Given that a particle has 
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a typical particular velocity against the liquid background of vo — which it acquires in 
collisions, see below — the viscous length ( is the typical distance that a particle has to 
travel before it again acquires the velocity of the background. The length ( may be estimated 
in the following way. From the equations of motion f<y = m,v, i.e., 

Vy = — ^T Vy > ^ 
we obtain a simple exponential decay of any excess residual y component of the velocity, 

67^77 



v y (t) = v exp 



-t 



rrii 



(5) 



with a time constant r = mj/67rrj?7. A typical excess velocity vq created by a collision of two 
particles 1 and 2 is v « + r 2 )7- The product 

^ v r = + r 2 /r 1 ) A f/6nr] (6) 

estimates the viscous length for particle type 1 after a collision with a type 2 particle. The 
value of d is largest for r 2 /ri = rb/r s and we call it simply ( in the subsequent text. A 
detailed discussion is presented in App. [A| 

We propose the following physical picture of the segregation process: The particles move 
on average with the velocity of the liquid at their centers. The collisions between particles 
tend to drive the system to a state of a random particle distribution ||. The ratio of the 
viscous length to a typical inter-particle distance is a measure for the efficiency of this process 

- the larger the viscous length, the more effectively will a collision disperse the state of 
the system. Conversely however, in a highly viscous environment with ( < 1, collisions will 
create order. 

To see the involved mechanism let us first look for stationary states of a monodisperse 
system, which we will define as those states not leading to particle displacements in y- 
direction of more than half a particle diameter. Each single particle defines, due to its 
finite extension, a horizontal "lane" in the system, i.e., the area that would be covered as 
time passes if no other particles were present in the system. A second particle will undergo 
collisions with the first particle if (i) their lanes overlap and (ii) if their ^-coordinates differ 

- under periodic boundary conditions the particles cannot escape each other. 

For a finite system with a given height one can find a maximum particle concentration 
below which all particles can be in disjoint lanes. But for an infinite system this concentration 
tends to zero, since collisions will necessarily occur if the particle arrangements are random. 
On the condition that the particle area fraction does not exceed the maximum packing 
fraction of disks in a strip one diameter wide, cq = 7r/4, we can imagine a stationary, 
collisionless configuration of particles. Let the particles be arranged one after another in 
a horizontal lane, all with the same y-coordinate. Then let these lanes be stacked in y- 
direction. Varying the horizontal distances between particles and the vertical distances 
between lanes any concentration c < cq can be achieved. 

Such a particle arrangement seems to be highly singular, but it can under certain con- 
ditions be stable against perturbations. Let the vertical distance between lanes be small 
with respect to r = 1 and imagine one particle in one of the lanes being slightly displaced 
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in +?/-direction. The particle will either undergo a collision with another particle in the 
lane above, which will then reduce its y-coordinate again, or it will collide with a particle 
in its own lane, increasing the original displacement. If the viscous length ( is large 1), 
this collision may suffice to displace both collision partners to another lane (if at the same 
time also the mean free path is large enough, see App. [TJ). In contrast, if ( is small (^C 1), 
then the relative motion of the two particles is more akin to a sliding on top of each other, 
displacing each particle by ~ 1/2 in +y and — y-direction respectively. In a sufficiently dense 
system subsequent collisions with particles can the neighboring lanes can then restore the 
original vertical positions of the particles. 

In the polydisperse case we might now expect that mixed lanes, containing both big and 
small particles, are stable. However, arguments can be made against this supposition. If ( 
is small in comparison to the big particle radius then small particles in lanes containing big 
particles will be ejected from the big particle lane by collisions. The reason is that the big 
massive collision partner is not significantly displaced by one collision, whereas the small 
partner is. If then the neighboring lane contains small particles, it is not hard to absorb 
the additional expelled particle. If, however, the neighboring lane contains big particles, 
alternating collision series of the small particle with big particles below and above will 
establish an additional small particle lane in between the big particle lanes. 

Starting from a random configuration at small (, local groups of big particles will tend 
to align and expel small particles — neighboring groups may then join due to the stresses 
generated by the collisions with the small particles accumulating outside these groups. Fi- 
nally, stable lanes form segregates into a of particles in disjoint lanes. A typical simulation 
sequence of this ordering process is displayed in Fig. [|. 

To arrive at a quantitative description of the segregation process, we define an "order 
parameter" in the following way. Since we have observed a strong stratification of the flow 
into horizontal layers, we define for each particle species the area fraction in a horizontal 
strip of width of the average particle radius. For a completely segregated system, we expect 
the area fraction of small particles a s to be large whenever the fraction of big particles a& is 
small. Consequently, the quantity 

S = ([a b - a s - (a b - a s )] 2 ) 1/2 (7) 

is small for a random mixture of the two species and assumes a large value when the system 
is stratified. The brackets denote the average over all examined horizontal slices. 

In Fig. |5| we show the time dependence of 5 for different values of ( for constant overall 
area fraction c = 0.6 and constant shear rate 7 = 0.01. We clearly see an initially fast size 
segregation process which becomes slower and slower and finally saturates to a value 5^ that 
depends on (. Due to the rather low value of ( and consequently a strong non-ergodicity of 
our systems, the sample-to-sample fluctuations are large and values of 0.2<5oo are typical in 
samples of size 100 x 100. 

The decreasing dependence of 5^ on C, which demonstrates the mixing or destabilizing 
effects of large (, is shown in Fig. || The figure shows data obtained for different fluid 
viscosities and shear rates 7, but constant overall area fraction. The scatter is rather large 
due to the abovementioned sample-to-sample fluctuations. At large (, the segregation does 
not increase significantly over the initial value. In fact, if ( is larger than the mean free path 
between particles, then the spatial distribution of particles does not differ much from that of 
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the corresponding, inelastic, sheared hard-core gas 0. However, at small £ the friction with 
the liquid is very large and causes ordering in bands oriented parallel to the flow containing 
alternatingly big and small particles. For a proposal to scale not only with respect to 7 and 
rj but also incorporating the concentration c, see App. |B|. 



IV. POISEUILLE FLOW 

To see the influence on the ordering effect of a different flow profile, we have exchanged 
the simple linear shear profile of Eq. ([|) by a quadratic Poiseuille profile, 

u(y) = 4L(| + i)(|-i)tw^. (8) 

We keep periodic boundary conditions in x-direction while introducing reflecting walls at 
the bottom and the top of the cell, where the flow velocity is zero in accord with no-slip 
boundary conditions. The shear rate now depends on the ^-coordinate in the flow, 

ly = 8j-v max . (9) 

Its modulus \ (d/dy)v x \ is largest at the walls and zero in the center of the flow. 

We have studied systems of concentration c = 0.6 and c = 0.4 with viscosity rj = 0.01 
and shear rates of 7^/4 = 0.01 at y = L/4 so that a typical zone in the flow corresponds to 
the crossover regime where ( ~ 1. After a long time £7^/4 ~ 100 we find that the particles 
concentrate in the central region of the cell, where the shear rate is smallest, cf. Fig. [7]. 
Close to the wall, and although the shear is largest there, we still find stratification — small 
particles gather at the wall. 

In a real experimental system, the liquid is strongly coupled to the particles. The momen- 
tum exchange between the two phases alters significantly the flow pattern, both on the scale 
of the particle radius, but at least in the viscous regime also on the scale of the container. We 
have tested our assumption of a laminar profile on the scale of the container by comparison 
to the results of a newly developed algorithm whose details are described in Refs. ||T0| , |IT 



There, we couple the particles to the liquid by Stokes-type drag forces proportional to the 
particle radius and the local velocity difference to the flow field — corresponding to Eq. (^|). 
The drag is introduced as a point-force term in the Navier-Stokes equations for the liquid 
surrounding the particles and accounts for the momentum transfer from the particles to the 
liquid. The Navier-Stokes equation is solved by a finite difference scheme, where the grid 
has a resolution on the order of the particle size. This technique is appropriate to observe 
the long range effects of the hydrodynamic interaction between the particles. The result of 
this simulation is shown in Fig. [8|, for the same choice of parameters as for the fixed velocity 
profile. 

Since there is no drastic visible difference to the fixed profile calculation shown in Fig. |7|(b) 
with the same physical parameters, we have a posteriori justified the choice a fixed liquid 
profile. To be more quantitative, we have checked the velocity profile of the combined 
particle-liquid system against the theoretical expectation for laminar liquid-only flow (see 
Fig. ^|) with constant viscosity. For moderate Reynolds numbers (~ 100) we find that there 
are only very small deviations from a parabolic flow profile. However, the particle load 
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increases the effective viscosity of the mixture and thus the overall flow velocity is reduced 
in comparison to the case of the pure liquid. In Fig. ^ we recognize a slight flattening of 
the parabolic profile near the center of the flow. Such a behavior is characteristic for shear 
thinning liquids and is here caused by the locally increased viscosity in the flow center, 
corresponding to the increased vertical momentum transport in the denser central region. 
At these Reynolds numbers, we measure only small velocity fluctuations of the liquid along 
the horizontal direction, i.e., the profile is stationary. 



V. CONCLUSION 



We have studied sheared bidisperse granular systems under conditions of simple shear and 
Poiseuille flow. We find that the presence of the liquid can induce particle size segregation. 
The criterion to decide whether this segregation mechanism will occur is that the ratio of 
shear rate to viscosity must be small. This statement is equivalent to saying that the viscous 
length in the system should be small when compared to a typical linear scale in the problem 
as, e.g., the particle radius. The particles arrange themselves in bands moving with the flow 
that contain alternatingly big and small particles. 

Also in Poiseuille flow we observe the formation of bands of particles of different size. 
Here, in addition, the variety of time scales in the flow leads to an aggregation of particles 
in the zones of low shear rate. These results have been verified against simulations using a 
full Navier-Stokes simulation for the liquid. 

The proposed segregation mechanism relies on the presence of a liquid phase. It is 
thus very different from known mechanisms in dry granular media, where gravity induced 
avalanches occur and separate particles species whose static angles of repose differ. 

It should be very interesting to study the behavior of the system with more than two 
particle species or a whole continuum of species or the dependence on the particle radius 
ratio. Moreover, three dimensional calculations are desirable and quantitative comparisons 
to experiments are necessary. 
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APPENDIX A: VISCOUS LENGTH 



We obtain an estimate (when no further collisions occur) for ( by integration of the 
equations of motion f^j = mjV. Taking the drag force from Eq. (^), we have obtained for 
the y component of the equation of motion (f|), which we print here once more: 

Vy = ^~ Vy - ( ) 

1 1 Li 
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This relation leads to a simple exponential relaxation of the initial excess velocity v y (0) 

67rriT) 



v y (t) = v y (0) exp 



-t 



rrii 



In analogous fashion, we obtain the expression for the x component, 

67rrj?7 r rt 



rrii 
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dt' Vy{t') 



(A2) 



(A3) 



which we can integrate by standard methods to find v x (t). We define the viscous length as 
the norm of the vector valued integral 



c 



dt [v(t) - u{y{t))\ 



(A4) 



and use the solutions of ( |A1|) and (|A3|) to obtain, 

v x (0) 



c 



rrii 



67rrj?7 



Vy(0) 



(A5) 



We then find typical values for vo by a consideration of a two particle collision — say 
between particles with label 1 and 2 — assuming that the initial velocities equal v, = u(?/j) 
according to their different y positions in the flow (for a more complete discussion of inelastic 
two particle collisions, see, e.g. [12]). The velocities after the collision, in the reference frame 
comoving with the liquid at the initial position of particle 1, are 



vi(0) 



rrii 



sin a cosa(l + e) 
sin 2 a — e cos 2 a + 1 



(A6) 



Here, b denotes the impact parameter and sin a = b/ (ri + r2), cf. Fig. |]. Thus, apart from 
order-one geometrical factors and some e dependence, the velocity of the scattered particle is 
13| m (/i/mj)(ri + r 2 )7. Therefore, ( becomes largest for the small particles after a collision 



involving a big and a small particle: 



/x(l +r b /r s )'j 

6711] 



(A7) 



APPENDIX B: AN ESTIMATE FOR THE MEAN FREE PATH IN THE 

MONODISPERSE SYSTEM 

The ratio of the viscous length to the mean free path of the particles in a non-viscous 
environment is probably an important dynamical characteristics of the system. If the mean 
free free path is much shorter than the viscous length, the additional background viscosity 
will not have a big effect. If on the other hand the mean free path is much longer than the 
viscous length, then the behavior of the system is viscosity dominated. Here, we would like 
to give an estimate of the mean free path in a monodisperse system for particles moving 
in the vertical direction. To allow for a simple calculation in the stationary situation, we 
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resort to a simple hypotheses for the system's configuration at large times: due to the initial 
disorder the systems arranges itself such that the number of horizontal particle lanes is 
maximum. 

We note that under these circumstances the average horizontal distance d between two 
particles is set by the packing fraction. If the number of lanes is maximum, their width must 
be the smallest possible, namely 2r. One particle covers an area of tit 2 within the available 
area 2rd. Consequently, the overall area fraction is 

7rr 2 1 r _ . 

C= = -7T-. Bl 

2rd 2 d K J 

As in the previous section, we now assume that we perturb the trajectory of one particle 
by giving it a vertical velocity of order rj, which is of the same order as the velocity difference 
between two lanes 2rj. If d = 2r, the system will not allow particles to penetrate into the 
neighboring lane. This situation is the densest packing compatible with a stationary state 
of the system, Co = 7r/4. As long as d < 4r, or equivalently, c>ci = 7r/8,a particle will 
only occasionally be able to pass a lane. Considerations of the particle geometry apart, the 
probability for a hit should be proportional to the time spend in the lane by the scattered 
particle divided by the average time between the pass of two successive particles in the 
neighbor lane. The mean free path i is given by the condition that this ratio be about 1, 
i.e., 

* ~ (B2) 



(d-2r)/2ir, 



or, 



ilr » (d-2r)/2r = -tic- 1 - 1 = — - 1. (B3) 
4 c 

For even lower concentration, the particle has a good chance to pass one or even several 
lanes, each with probability 1 — Pha ~ 1 — [2r/^r\/[(d — 2r)/27r] = 1 — {2r)/(d — 2r). The 
probability to survive a distance x/r without hits is hence distributed exponentially, 

p(x/r)~(l-j> hi t) a,/2r , (B4) 

which yields — by normalization and determination of the expectation value - 

£/r = -l/ki{l-p m ). (B5) 

For small concentrations (and thus also small hitting probabilities) this equation may be 
expanded to yield the same form as (TB3I) 



e/r ~ 1/pMt ~ ^ - 1. (B6) 
c 

It is interesting to note that i does not depend on the shear rate but only on geometrical 
properties of the system. This observation is in favor of our suggestion that the ratio of 
viscous length to particle radius £/r, as proposed in the main text and here in dimensional 
form, collapses the simulation data for 5 at a fixed given c. If 7 <C 1, we presume that Q/C. 
may be a good scaling variable, even at different area fractions. This may be an interesting 
question to investigate. 
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FIGURES 



FIG. 1. Sketch of the model geometry for shear flow. The simulation cell of size L x L is 
repeated in x-direction (periodic boundaries). The cell images in ±y-direction are shifted by an 
amount of ±jLt/2 to reflect the particle displacement at the top and bottom of the cell which 
grows linearly in time due to the shear (Lees-Edward boundary condition). 

FIG. 2. Geometry of the collision between two disks in the center of mass frame. The particles 
are assumed to initially move parallel to the flow, such that laboratory from and center of mass 
frame are aligned. 

FIG. 3. Simulation snapshots at non-dimensional time jt = (a) and 195 (b). The shear rate 
in this system is j = 0.01, the viscosity 77 = 0.001. 

FIG. 4. Simulation snapshots at non-dimensional time jt = (a), 20 (b), 70 (c), 110 (d) 
and 220 (e). The shear rate in this system is 7 = 0.01, the viscosity rj = 0.01, ten times larger 
than in the preceding figure. Different shades of grey indicate the modulus of the x- velocity of the 
particles. 

FIG. 5. Time dependence of the segregation parameter 5 for simulations with different viscosity 
r] = 0.001, 0.004, 0.01, 0.03 (bottom to top) corresponding to ( ~ 2.5, 0.63, 0.25, 0.08 [according 
to Eq. (|A7|)] (bottom curve to top curve) vs. dimensionless time jt on the abscissa. 

FIG. 6. Final values of the segregation parameter 5 plotted vs. viscous length £ [according to 
Eq. ( [A7D ] on the abscissa for several values of viscosity and shear rate in the system. Crosses denote 
computations with constant viscosity rj = 0.01, diamonds indicate constant shear rate 7 = 0.01. 

FIG. 7. Stationary particle configuration with imposed parabolic Poiseuille flow profile. The 
system size is 100 x 100, the particle volume fraction c = 0.4 (a), and 0.6 (6). 

FIG. 8. Final particle configuration in Poiseuille flow with velocity profile obtained by solution 
of the full Navier-Stokes equation with a point force term modeling the momentum exchange 
between fluid and particle phase. There is no difference in the physical parameters of the two 
systems. The system size is 80 x 40. 

FIG. 9. The x-component of the velocity profile along a cut in y-direction in the Navier-Stokes 
Poiseuille flow simulation at Re = 100. The solid line is the prediction for laminar flow without 
particles, the dashed line is a parabolic fit to the data points. 
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